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Abstract 

We investigate the behavior of the shear viscosity n{p) and the mass-dependent diffusion co- 
efficient D(m,p) in the context of a simple model that, as the crosslink density p is increased, 
undergoes a continuous transition from a fluid to a gel. The shear viscosity diverges at the gel 
point according to n(p) ~ (p c — p)~ s with s ~ 0.65. The diffusion constant shows a remarkable 
dependence on the mass of the clusters: D(m,p) ~ m -0 69 , not only at p c but well into the liquid 
phase. We also find that the Stokes-Einstein relation Dn oc fceT breaks down already quite far 
from the gel point. 
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I. INTRODUCTION 



When a system of polyfunctional molecules is crosslinked, the transport properties such 
as the shear viscosity and the diffusivity can be dramatically affected. In particular, the 
diffusivity decreases as the number of crosslinks is increased and the shear modulus increases, 
diverging at the critical point at which a gel is formed. The diffusivity remains finite as the 
system gels since monomers and small clusters can diffuse through the tenuous structure 
that characterizes the amorphous solid close to the critical point. Although gels have been 
studied for many years [lj, their critical behavior remains poorly understood. In particular 
the question of whether or not there exist universality classes into which different materials 
can be grouped remains largely unanswered. 

In this article, we report on extensive molecular dynamics simulations of a simple model 
for a gel. We study the system on the fluid side of the gel point from the simple liquid limit 
into the critical region. We investigate the structural properties of clusters and calculate 
both the shear viscosity r)(p) and the mass-dependent diffusion constant D(m,p) as function 
of the crosslink density p. We find that as p — > p c , rj(p) ~ (p c — p)~ s with s ~ 0.65, a value 
somewhat smaller than that conjectured by de Gennes Q] on the basis of an analogy with 
a random superconductor network and also predicted recently by Broderix et. al. [3| for a 
Rouse-like model network. The mass-dependent diffusion constant D(m,p) ~ m -°- 69 f or a 
range of p near the critical point and 3 < m < 50. This behavior is consistent with earlier 
results for p = p c jj] and rather close to a prediction made on the basis of a simple scaling 
argument. On the other hand, our value for s is somewhat lower than the one found by 
Kiintzel et. al. in a recent article in which s = 0.8 is found by a theoretical analysis 
of the Zimm model. The diffusion coefficient D(m,p) — > const, as p — > p c for m at least as 
large as 10 but displays critical behavior in the next leading term. It is also worth noting 
that, in contrast to simple liquids, the product D(p)r](p) is not a constant but rather reflects 
the divergence of r\ at the gel point. 

The structure of this article is as follows. In section |H] we describe our model and the 
computational details. Section lllll contains a discussion of the geometric properties of the 
clusters and the nature of the percolation transition. The shear viscosity calculation and 
results are described in section HVl and results for the diffusion constants are found in0 We 
conclude with a short discussion in section IVII 
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II. MODEL 

The model is similar to the one employed in [4], but we include the details below for 
completeness. Our system is composed of N = L 3 (L = 10, 13, 15, 20 and 30) particles 
interacting pairwise through the shifted Lennard- Jones potential 

U{v)=l (1) 
I otherwise, 

where Ulj(t) = Ae((a/r) 12 — (cr/r) 6 ). All of our simulations are 3D constant energy molecu- 
lar dynamics (MD) simulations corresponding to an average temperature of ksT/e ~ 1 and 
density $ = 0.8a~ 3 . These choices ensure that the system is in the liquid-phase region of 
the phase-diagram (f| ?]]. We use periodic boundary conditions and a time step of magni- 
tude dt = 0.005r, where r = ^Jmo 2 Je is the reduced Lennard- Jones time. From a typical 
equilibrium state of this liquid we let the particles form a specified number n of permanent 
chemical bonds if they come closer than .12, coinciding with the minimum 

of U(r). The bond interaction is a harmonic oscillator potential t/harm( r ) = 1/2/cr 2 : in our 
simulations we take ka 2 /e = 2.0 (different from ^J). Note that this way of adding bonds 
violates energy conservation; indeed we actually pump energy into the system when adding 
bonds. To compensate we cool down the system again after having established the required 
number of bonds. With this bonding procedure cross linking is very fast — the average dis- 
tance between the particles is comparable to r c , so a large number of particles are available 
for bonding at any given instant. Each particle can bond to a maximum of / = 6 other 
particles (excluding itself), and the cross link density p is then given in terms of the number 
of bonds n as p = 2n/ f'N. Any number of particles, if fulfilling the conditions above, can be 
cross linked per time step, but we halt the bond formation when p reaches a predetermined 
value. 



III. GEOMETRIC PROPERTIES 



Before discussing the dynamic properties of this model, we need some basic information 
about the static properties. In this section we determine the geometrical percolation point p c 
as well as the two critical exponents u, the correlation length exponent, and 7 the exponent 
characterising the divergence of the weight average cluster mass (of finite clusters). We 
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FIG. 1: Fraction of systems W(L,p) percolating in the x-direction as a function of p and for five 
different system sizes as indicated on the plot. The lines are guides for the eye, except in the case 
L = 30 for which the data are fitted to a stretched exponential [llj]. We estimate p c = 0.2565. 

follow a similar procedure to the one used and outlined in DDD . In order to find p c we 
calculate numerically the fraction W(L,p) of percolating systems of size L 3 with a bond 
density p. This function is plotted in Fig. ^ for all five system sizes. The crossing points 
of the different curves seem to coincide, and the corresponding value of p is thus a good 
estimate of p c jj, Q|: From the figure we determine p c = 0.2565 as in Finite size 
scaling theory predicts that W(L,p) does not depend on L and p separately but only on the 
combination L/« (and the sign of o- Pc ) where { = |p - ^~ is the correlation length and 
v the correlation length exponent [12J. Thus we may write 

W(L,p) = f(L^(p-p c )), 



(2) 



where f(x) is a scaling function. To test this hypothesis we replot the data for W(L,p) from 
Fig. ^ in Fig. |21 as a function of L l l u (p — p c ) with p c = 0.2565 as determined above and 
v = 0.9. The collapse is very good confirming the correctness of the values for p c and v. 

To compare with percolation theory we need one more exponent, and here we consider the 
behavior of the weight average cluster mass M w . In the thermodynamic limit the expected 
behavior is M w {p) ~ \p— p c \ ~ 7 [12]: Therefore we compute M w as a function of p for different 
system sizes, and in Fig. El we plot the results in the form M w / L 1 ^ versus L l l u (p c — p) with 
7 = 1.8 being the expected 3D percolation value and v and p c as determined previously. 
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FIG. 2: Same as Fig.^ except here W(L,p) is plotted as a function of L 1 ^ '(p—p c ) with p c = 0.2565 
and v = 0.9. The data collapse very nicely in agreement with finite size scaling theory. 
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FIG. 3: Scaling plot of the weight average molecular weight M w . The quality of the data collapse 
confirms 7 = 1.8 in accordance with the 3D percolation value. 

Again there is a very nice data collapse, and we therefore conclude that here as in 0,0,0] our 
system is consistent with the 3D percolation universality class in so far as static properties 
are concerned. 
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IV. VISCOSITY 

We measure the shear viscosity rj(p) by using the appropriate Green-Kubo formula 

Q: 

1 f 00 

ri = y^J^ dt(a xy (t)a xy (0)), (3) 
where V is the volume and <y xy {t) the xy component of the stress tensor: 

N N 

a xy{t) = ^2 mV *,i V y,i + ~ Vi)fxM- ( 4 ) 

i=l i=l j<i 

In this equation f x ^j is the x component of the force from particle j on particle i, and 
the meaning of the remaining terms is self explanatory. In the simulations we average over 
several hundred samples for each p and over three off-diagonal components (xy, yz and zx) 
of the stress tensor to obtain slightly better statistics. It is important to note that we have 
discarded any sample containing a spanning cluster since for such a system the viscosity is 
not defined, i.e. the right hand side of Eq. @ diverges. Although we have simulated very 
long runs (up to t = 750r) the stress correlator C aa (t) = (vxyi^Cxyify) has still not decayed 
completely, and it is necessary to add by hand an additional contribution, in particular for 

< c < 0.3 seems to fit 
to believe that this 
is the appropriate form. See |9| for a thorough discussion of this point. 

In Fig. |3] we have plotted the resulting values for the viscosity for different systems sizes 
and at different stages of the cross linking. We note the clear power law behavior outside 
the critical region, and a fit to the L — 10 data in this region yields s = 0.65. The line 
7] oc (p c — p)~ 0M5 has also been drawn on the plot, and it is apparent that the data are 
consistent with this exponent. For large p, p > 0.23, there are larger error bars and this 
will also affect the scaling plot. Since the viscosity diverges at the critical point with an 
exponent s > 0, the finite size scaling form is 

V (p,L) = L- s ^g(Z/L) p<p c , (5) 

where g is a scaling function with the limits 

x -s/u x _^ o 

g(x) oc { (6) 
const, x — > oo, 



p close to p c . A stretched exponential C aa (t) = aexp(—bt c ) with 0.] 
the data well for lone times, and there are also theoretical reasons 15] 
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FIG. 4: The dimensionless shear viscosity as a function of p c — p for different L. 
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FIG. 5: Same as Fig. but here plotted in a scaling form with s = 0.65. 

and £ ~ (p c — p) - ^ is the correlation length, c.f. Sec. IIHI Therefore we plot in Fig. \^i]L s ^ u 
versus L 1 '"(p c —p), and the collapse is quite good outside the critical region with s = 0.65, 
whereas there is a larger scattering of the points for p closer to p c . 

de Gennes has suggested Q a value of s ~ 0.7 based on an analogy between gelation 
and conductance in a random mixture of normal and superconducting elements, and nice 
agreement with this was found in a related model in Here we have observed a slightly 
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smaller value for s. 



V. DIFFUSION 

n 

In this section we extend our earlier study 4] on the diffusion of clusters. Previously we 
were concerned mainly with the behaviour of the diffusion constant D(m,p) as a function 
of cluster mass m at the gelation point p = p c . Here we address the p dependence of D 
for different clusters and the validity of the Stokes-Einstein relation D(p) oc ksT /f](p) for a 
given cluster mass. We restrict our attention to the L = 20 system. 

To determine the diffusion constant we use the Einstein relation: 

±((r m (t)-r m (0)T}^D(m,p), (7) 

where r m (t) is the center-of-mass position of a cluster of mass m at time t, for a given value 
of p (for clarity of the presentation we omit the explicit dependence on p in the notation). 
When calculating the diffusion constant numerically we have averaged over all clusters of a 
given mass m and over several hundred crosslinkings, and we have discarded any percolating 
samples. This has been done mainly for consistency when comparing with 77, but in any 
event we do not expect this to affect the diffusion of any but the very largest clusters. 

First we examine the convergence of Eq. (J2J) by plotting in Fig.Elthe behavior of ((r m (t) — 
r m (0)) 2 )/6t for m — 1 (monomers) as a function of time and for three different values of 
p. From these curves we clearly see the existence of long-time tails in the velocity auto- 
correlation function. Consider the "Green-Kubo" formula corresponding to Eq. (|7jl: 



((r m (t)-r m (0)) 2 ) 



/ ds(v TO (s).v m (0))(l-s/*)- ( 8 ) 
Jo 



t 

The dominant contribution to ((r m (t) — r m (0)) 2 )/t at large times is 

((r m (t) -r m (0)) 2 ) =D(TO)p) _ HdsC^is). (9) 



t 

where Cw (s) = (v m (s) • v m (0)) is the velocity auto-correlator. Therefore, a power law tail 
ci™\s) ~ t~ a in the velocity auto-correlation function will translate into a corresponding 
power law tail ((r m (t) — r m (0)) 2 )/£ ~ D(m,p) + const. t l ~ a in the Einstein relation. In 
simple liquids a value of a = 3/2 is ubiquitous [1J|, and has also been observed for gelating 
systems in [{J: here we find that the same power-law provides a very good fit to the data 
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FIG. 6: ((?"i(t) — ri(0)) 2 )/6t a function of time for monomers, and for three different values of p: 
0.125, 0.2 and 0.25 from top to bottom. The long time tails are clearly visible, and the solid lines 
are fits to the same functional form (see text). 
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FIG. 7: Same as Fig. H3 but this time for clusters of size 10. 



for all m, but in particular for the small clusters. In Fig. El we have also plotted these fits 
to the power-law a + bt^ 1 ^ 2 , and the deviation from the simulation results at early times is 
barely visible. In figure Fig. [7| we have done the same for clusters of mass 10, and we see 
the same behavior. The agreement is slightly worse, presumably due to poorer statistics of 
larger clusters. We note the existence of a maximum in all of the curves (though not visible 
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on Fig. IHlfor m = 1) for ((r m (t) — r m (0)) 2 )/t. By differentiating Eq. (jHJ) this can be shown 
to occur at t m , where t m is the solution to 

p t>m 



dsC^\s)s = 0. (10) 
An obvious consequence of the fact that for t > t m (using Eq. (jHJ)) 

| ((r - (t) - r - (0))2) ^^)(t)<o (11) 

is that Ct\t) becomes negative (anti-correlation) for large t and stays negative thereafter. 
This means that the l/\/t tails in Figs. EH3 correspond to a negative t -3//2 tail in Cm (t). 
We also note that t m is an increasing function of m and a decreasing function of p. 

The error made by taking D(m,p) to be the value of {{r m {t) — r m (0)) 2 )/6t at the end 
of the simulation time t = 120r is neglibly small: for p = 0.2 we have compared with 
simulations that are twice as long, and at least for the 20 lightest clusters for which we had 
good enough statistics, the error was less than 5%. For the smallest m where the statistics 
are very good, one can also obtain D(m,p) from a power-law fit as mentioned above, and 
the outcome is still consistent with the statement just made (the error here is even much 
smaller than 5%). 

In jj] we studied D(m,p c ) and we found the power-law D(m,p c ) ~ m~ om . We have 

repeated this study up to clusters of size 50, and we observe the same behavior over the 

entire range. For p < p c , we see the same power law as a function of m, at least for small 

cluster sizes. The quality of the statistics for larger cluster sizes is insufficient to determine 

whether there is a cross-over or cut-off as m — > m*(p), where m*(p) ~ (p c — p)~ 1 l a is the 

mass of the largest cluster, but it seems likely that there is. de Gennes has argued that for 

masses 1 < m < m*(»), Dim) ~ m -( u + s )/(f 3 + / y) Q n the basis of a Stokes-Einstein relation 

□ 

with a mass dependent viscosity |2J. Here f3 is the exponent that describes the decrease of 
the order parameter near percolation: x ge \ ~ (jp — PcY, where x ge i the fraction of particles on 
the spanning cluster and p — > p c +- The other exponents have been introduced already. By 
using the appropriate scaling relations for 3D percolation 12], the exponent can be rewritten 
so the prediction is D(m) ~ m~ 2 ^ u+s ^^ du+1 \ where d = 3 is the Euclidian dimension. With 
our values for the remaining exponents we get: 

D(m,p) ~ m~ - 69 for 1 < m < m*(p). (12) 
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FIG. 8: Diffusion constant as a function of p for clusters of three different sizes: m = 1, m = 2 
and m = 3 from top to bottom. The solid line is a fit to the function D(p) = a(p c — p) b + D c and 
D c = 0.0398, a = 0.131 and b = 1.103 (see text). 

This is in very good agreement with our simulation results within the observed power law 
regime. The theoretical prediction can be rewritten as D(R g ) ~ Rg 1+S ' u ^ where R g ~ m l l D s 
is the radius of gyration and Dt the fractal dimension. This form of the relation has 

n 

sometimes (see for example |16() been used to infer s from the scaling of D with R g , but to 
the best of our knowledge the present study presents the first direct verification of such a 
link. 

However, even in this regime one would expect some additional p dependence of the 
diffusion coefficient, a point not addressed in To this end, we plot in Fig. \^D(m,p) as 
a function of p for monomers, dimers and trimers, and we see that the diffusion constants 
decrease (almost linearly) as a function of p. Moreover the curves seem to fit nicely to the 
functional form D(p) = a(p c — p) h + D c , with a value of the exponent b = 1.1. In Fig. El we 
have made a similar plot for masses m = 2 ... 10, and the trends observed above appear to 
carry over to larger masses. The curves are roughly parallel, and therefore it is not unlikely 
that the value of b is independent of m, but we are unable to confirm this from a fit to the 
data: the exact value of the exponent appears to be very sensitive to noise in the data. 

Finally we demonstrate a striking violation of the Stokes-Einstein relation when ap- 
proaching the gelation transition. The idea that D oc 1/rj is used so widely that one may 
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FIG. 9: Same as Fig. El but for clusters of sizes m = 2 ... 10 from top to bottom. The solid line 
is a fit, and here D c = 0.0064, a = 0.0324 and b = 1.029. 
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FIG. 10: Plot of D(m,p) times rj: obviously this is only approximately a constant, as predicted 
by the Stokes-Einstein relation, far away from p c . 

sometimes forget its lack of universal validity. In Fig. fTUl however it is clear that D(m,p)r) 
increases significantly when p p c . This is consistent with our previous observations that 
whereas r\ diverges at the gelation point, D(m,p) approaches a non- vanishing constant even 
for large masses m. Further away from the gelation point p < 0.20 there does however seem 
to be an approximate proportionality between D(m,p) and r\. However, in Figs. IHHHlwe saw 
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indications that D(p) ~ a(p c — p) b + D c with b > 1 whereas 77 ~ (p c — p) 065 , and so this 
apparent proportionality is at best only approximate. 

VI. CONCLUSIONS 

The main results of this study are the power law behavior of the mass dependent diffusion 
coefficient which seems to hold well away from the critical point, the failure of the Stokes- 
Einstein relation and the result s ~ 0.65 for the critical exponent of the shear viscosity. This 
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seems to support the conjecture 



last result, taken together with other recent results 
that the gelation transition is not classifiable in terms of a single universality class: exponents 
in the range 0.3 < s < 0.7 have been found for models that seem, on the surface, to be very 
similar. The experimental situation also does not provide much evidence for universality: 

both exponents near s = 0.7 17|, [l8|, [l9| and in the range 1.1 < s < 1.3 0, 0,0,0,3 



have been reported. We sound a note of caution here: The determination of exponents 
through finite size scaling is not very precise, especially when quantities that are as difficult 
to calculate as the shear viscosity form the data set. However, it seems very unlikely that 
the errors are large enough that a factor of more than 2 in the exponent could be explained 
that way. 

The mass-dependent diffusion coefficient in this model displays a power law behavior 
D(m,p) ~ m -°- 69 consistent with a scaling argument of de Gennes j2|. Reexpressing this 

Dt 

in terms of the radius of gyration of clusters through m = R g where Df = 2.5 is the 
fractal dimension of the percolating cluster, the scaling prediction is D(R g ,p) ~ R g ^ 1+S ^ w \ 
This yields an estimate s = 0.65 for the viscosity exponent, in good agreement with the 
direct calculation from the Green-Kubo formula. Whether this connection between diffusion 
and viscosity is general or specific to the present model and whether there exist a similar 
relationship between diffusion and the elastic shear modulus in the solid phase remains a 



subject for further study. Using a quite different model, del Gado et al. [16J] have studied 
the self diffusion of crosslinked polymer clusters on a lattice by bond fluctuation dynamics. 
They have also used this scaling ansatz to infer the critical exponent of the shear viscosity 
and found s ~ 1.3. Their result translates to a mass dependence of the diffusion constant 
D(m,p) ~ m -1 , very different from that of the present model. 

Finally, we have shown that as the fluid becomes more viscous there is a breakdown of the 
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Stokes-Einstein law Drj oc k B T that generally holds for simple liquids. For relatively small 
concentrations of crosslinks, this product varies only very little but for p fa p c the divergence 
of the viscosity begins to dominate the diffusion constant which seems to saturate for all 
cluster sizes studied at p c . 
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